/•""\ 



MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
ARTIFICIAL INTELLIGENCE LABORATORY 



A.I. Memo No. 440 September 1977 



DENSITY RECONSTRUCTION USING ARBITRARY RAY SAMPLING SCHEMES 



Berthold K. P. Horn 



ABSTRACT, 

Methods for calculating the distribution of absorption densities in a 
cross section through an object from density integrals along rays in the 
plane of the cross section are well-known, but are restricted to particular 
geometries of data collection. So-called convolutional-backprojection- 
summation methods, used now for parallel ray data, have recently been ex- 
tended to special cases of the fan-beam reconstruction problem by the addi- 
tion of pre- and post-multiplication steps. In this paper, I present a 
technique for deriving reconstruction algorithms for arbitrary ray-sampling 
schemes: the resulting algorithms entail the use of a general linear opera- 
tor, but require little more computation than the convolutional methods, 
which represent special cases. 

The key to the derivation is the observation that the contribution of 
a particular ray sum to a particular point in the reconstruction essentially 
depends on the negative inverse square of the perpendicular distance from 
the point to the ray and that this contribution has to be weighted by the 
ray-sampling density. The remaining task is the efficient arrangement of 
this computation, so that the contribution of each ray sum to each point 
in the reconstruction does not have to be calculated explicitly. 

The exposition of the new method is informal in order to facilitate 
the application of this technique to various scanning geometries. The 
frequency domain is not used, since it is inappropriate for the space- 
variant operators encountered in the general case. The technique is il- 
lustrated by the derivation of an algorithm for parallel-ray sampling with 
uneven spacing between rays and uneven spacing between projection angles. 
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BACK6R0UND AND MOTIVATION. 

Recent interest in computerized axial tomography as a means of deter- 
mining absorption densities in a cross section through an object has led 
to a variety of basic algorithms [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Part 
of this interest stems from the diagnostic benefits derived by the medical 
community from scanners utilizing X-ray sources which provide cross sections 
of the head, body and, soon, the heart. Reconstruction from the mass of 
data generated by many ray samplings was not feasible before the advent of 
small, fast computers, and the choice of reconstruction method depends to a 
large degree on the speed with which such a computer can perform the calcu- 
lations. As a result the so-called convolutional-backprojection-summation 
algorithm has emerged as the method of choice and is largely displacing 
competing methods using two-dimensional Fourier transforms or iterative 
solution techniques used to solve large sets of sparse equations. These 
other methods do still find application in specialized areas where speed 
is not the main criterion of success. A further advantage of the methods 
based on convolution is that each collection of density integrals, also 
called a projection, can be treated in a separate computation [6, 11], 

Reconstruction methods developed so far, however, have mostly been 
suited to the parallel -ray projection method of data collection, commonly 
employed in early, slow computerized axial tomographic scanners [9], Here 
density integrals or ray sums are sampled evenly along a line perpendicular 
to the rays (see figure 1); such a collection of data is called a projection, 
and projections are formed for a set of projection angles evenly spaced 
over either 180° or 360°. Reviews of a variety of reconstruction algorithms 
^ for this ray-sampling scheme may be found in several references [12, 13, 14]. 
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Since X-rays cannot be focused or deflected as visible light rays can, 
the pencil beams used for parallel -ray sampling are obtained by tight 
collimation of radiation emitted from an X-ray source radiating into a large 
solid angle. Most of the output of the source is therefore wasted. Since 
a certain number of X-ray photons must be absorbed in order to get a suffi- 
ciently accurate estimate of the density integral along the ray, a great 
deal of time elapses before all ray sums have been observed. In the mean- 
time, the object may have moved. For these and other reasons, modern 
scanners use fans of rays striking a multiplicity of detectors (see figure 
2). A whole projection may now be measured in the time it would have taken 
to measure a single ray sum with the older system [15, 16]. 

One difficulty with the so-called fan beam approach is that ray sums are 
no longer evenly spaced in terms of ray direction and distance of rays from 
the center of the region being scanned. As a result, conventional reconstruc- 
tion techniques do not apply without modification. Resorting the ray sums 
and interpolating to approximate parallel-ray data has not proved very ef- 
fective, because accuracy is compromised by the interpolation step [15, 16, 
17, 18]. 

Convolutional reconstruction methods have been modified, however, to 
deal with two very special cases of this ray-sum collection scheme [19, 20, 
21]. The first method applies to the situation where the fan is sampled 
evenly along a line at right angles to the line connecting the source to 
the center of the region being scanned. Such data collection can be achieved 
only with a detector array that co-rotates with the source of radiation. 
This puts a demand for exceptional stability on the central detectors in 
the array, since points near the center of the region being scanned are 
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"seen" only by a few detectors during the complete scan [22]. 

The second method applies to the situation where the detector array 
lies on a circle about the center of the region being scanned with radius 
equal to the radius of the circle on which the source moves. This geometry 
lends itself to the use of a fixed detector array with consequent simplifi- 
cation of the scanner mechanics. Since they lie on the same circle, there 
is a spatial conflict between the source and the detectors. If they are 
placed on circles with differing radii, the special case solution no longer 
applies. The latter geometry is in fact common amongst proposed fast scanners, 

Clearly a method is needed for deriving algorithms similar to these 
modified convolutional methods for data collected by arbitrary sampling of 
the ray-sum space. Unfortunately, as it turns out, convolutional-backpro- 
jection-summation techniques apply only to a few special geometries. Even 
/**\ the two fan-beam reconstruction methods mentioned above augment the convolu- 

tional step of the algorithm with a premulti plication of each ray sum by a 
factor depending on the position of the corresponding ray in the fan. Further- 
more, both involve the use of a postmulti plication during the summation step 
with a factor which depends on the position of the point being reconstructed 
relative to the fan currently being treated. 

While the main impetus for this work comes from the computerized X-ray 
transverse axial tomography application, similar methods are of importance 
in such other fields as radio astronomy [2, 3] and electron microscopy [4, 5]. 
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PREVIEW. 

The algorithms developed here use general linear operators. Operations 
using general linear operations can be thought of as spatially varying con- 
volutions, where the "kernel" or "point-spread-function" is allowed to de- 
pend on the position at which the operator is applied. The derivation de- 
pends on the following observations, which will be elucidated in the next 
few sections: 

«3f The contribution made by a particular density 
integral or ray sum to a particular point in the 
reconstruction is a function of the perpendicular 
distance from the point to the ray. 

•X* This contribution is essentially proportional to the 
negative inverse of the square of the distance, ex- 
cept for rays passing very near to the point in 
question. 

«&The contribution of a particular ray sum has to be 
divided by the local ray-sampling density, to account 
for uneven sampling of the ray-sum space. 

•JfThe ray-sampling density is simply the inverse of the 
Jacobian of the transformation from a convenient uni- 
form scanning coordinated system to the coordinates 
used in parallel ray reconstruction. 
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•X* Using a general linear operator, it is possible 
to arrange the computation efficiently for most 
scanning geometries of interest. That is, each 
generalized projection gives rise to a separate 
computation and it is not necessary to determine 
the contribution of each ray sum to each picture 
cell explicitly. 

«$£ For a few special cases, the general linear operator 
is spatially invariant and thus is simply a con- 
volution. Parallel -ray sampling is the best known 
example of this. 
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SOME PRELIMINARY DEFINITIONS . 

The notation used here is similar to that used by Lakshminarayanan 
[19, 6]. The set of rays sampled is a finite subset of the two-parameter 
family of straight lines in the plane. Various ways can be envisioned for 
designating particular rays. We may, for example, specify the inclination 
e of a ray (relative to the upright axis in figure 3), as well as the per- 
pendicular distance % from the center of the region being scanned. For some 
scanning geometries, other parameters will be more suitable, but for para- 
llel-ray systems this method is convenient, because, for this case, the 
projections correspond to evenly spaced values of e, while rays within a 
projection correspond to evenly spaced values of i. 

Let p(£,e) be the density integral or ray sum along the ray U,e). In 
practice we will be given only a finite set of these density integrals, 
corresponding to discrete values of i and e which depend on the scanning 
geometry. If we choose to use polar coordinates (r,<f>) to designate points 
in the region scanned, and let f(r,<|>} be the absorbing density at the point 
(r,<j>), then our task will be to reconstruct values of f(r,<j>), given a set 
of values of pU,e). 

One important quantity we will need is the perpendicular distance, t, 
from a given point to a ray. Using figure 3 again, we get, 

t = i - r cos(e - <j>) (1) 

If we let V be the value of % corresponding to t = 0, the case of a ray 
passing directly through the point, then i x = r cos(e - $), and so 

t = l - £' (2) 
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RADON'S FORMULA RE VISITED, 

The earliest known solution to the reconstruction problem is given by 
Radon in his paper of 1917 [1]. His result will not be rederived here, since 
advanced mathematical concepts are needed and because he has given such a 
clear account of the proof. To apply his formula, we have to assume that 
f(r,<j)) is bounded, continuous and zero outside the region scanned. Then 
pU,e) will also be zero outside a certain range for i. Further, pU,e) 
will be continuous. Now assume that the partial derivative of p(£,e) with 
respect to i is continuous, too. Radon's inversion formula then is [21, 15, 1], 



f(r,o) --7/ f V \) ^PU.e) d * ^ (3) 



-^ The above result does not strictly apply if some of the conditions — 

particularly the one regarding the continuity of the partial derivative of 
pU,e) --are violated. We may expect certain artifacts or reconstruction 
errors in and near regions where the assumptions fail to apply. The mag- 
nitude of the resulting errors depends on the details of the numerical 
approximations made to the above equations. 

The inner integral is singular, since t = 0, when 1 = i l (equation 2). 
This singular integral may be interpreted as 

£'-£ +00 

^/ (" i> & P<M> -* + "" f <-i>£p(i.e>d« (4) 

-°° £ £'+£ 

Integrating both terms by parts, we get 
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lim 
e->0 



1 pU'-e.e) + 1 pU'+e.e) + /(- \)pU,e) ds. 

|t| > e 



(5) 



Since p(£,e) is assumed to be continuous with respect to a, we can rewrite 
this 



Tim 

£+0 



/ F g (t) p(£,6) d£ 



(6) 



where 



F (t) -^ 



for |t| < e 
for |t| > e 



(7a) 
(7b) 



F^^\ 



Combining the above results, 



Tim 



-, r 11m f 

f(r,e) = -i- / I f (t) p(A,e) d£ de 

4lT 2 ^ n ^ 



(8) 



% Clearly each density integral or ray sum pU,e) 
contributes to each point in the reconstruction 
according to its distance from that point. 

5jc In fact, all but those rays passing wery close to 
the point, contribute with weight proportional to 
the negative inverse of the square of the distance 
from the point. 

^c The weight of the contributions of rays passing 
near the point is such that the sum of all weights 
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is zero. That is, 
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F(t) dt = (9) 



e 



The above three observations are important in the derivation of the new 
algorithms. The only other problem that will have to be tackled concerns 
the calculation of ray-sampling density. Then the techniques developed here 
can be applied to particular scanning geometries. 
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REVIEW OF THE CONVOLUTIONAL-BACKPROJECTION-SUMMATION METHOD. 



Let the inner integral (equation 6) discussed in the previous section 
be called gU.'.e). Clearly it can be thought of as a convolution of the 
original projection data pU.e) with the filter function F (t), since 
t = £' - i and F (t) are symmetric: 



e 

lim f 

gU 1 ,e) = e ^ J F £ (z l - i) pU,e) it 



(10) 
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We can then use the outer integral (equation 8) to calculate the densities 
from this convolved or filtered data: 

/2tt 
gU'.e) de (11) 



In practice we know only a finite number of ray sums and consequently have 
to approximate both of the above integrals by finite sums. If we choose to 
observe M projections evenly spaced in angle from e = to e = 2tt, we may 
approximate the outer integral (equation 11) by 

M - 1 
f(r,8) - -L- V g (£-) 69 (12) 

47F j = 

where 59 = (2tt)/M is the angular increment between successive projections. 
Note that QjU 1 ) is the convolved projection data of the j projection evalu- 
ated at V = r cos(e-<j>). For reasons of computational efficiency, we cal- 
culate the convolved data at only a small number of places — typically the 
same ones for which projection data is available. This allows the use of a 
single convolution per projection, independent of the location of the points 
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in the reconstruction. The value of g.(A'), needed in the above summation, 
must however then be estimated by interpolation from the values at those 
places where the convolution wa£ actually computed. 

This convolution will be discussed next. If we let W be the width or 
diameter of the region being scanned, and N the number of evenly spaced rays 
across this width (sampled by the detectors) we can approximate the inner 
integral (equation 10) by 



N - 1 

1 = o 



9 i'j = lL F i'-i p ii " (13) 



where si - W/(N - 1) is the uniform interval between successive rays in a 
projection and p... is the ray sum for the i th ray in the j th projection (see 
figure 4). Now the i' ray passes at a distance £' = i' S£ - w/2 from the 
origin, so this is the value of a' associated with g., .. From this relation- 
ship one can determine which values of g. ,. should be used in the interpola- 
tion for estimating g.(i'). One uses g . , . and g,., + }) . where 

i' = [(£' + W/2)/6iJ ( 14 ) 

We next turn our attention to the discrete approximation to F (t) 
(equation 7), 

w k 



o " E F k ( 15b ) 



F_ = -2 

k = 1 



The value of F Q is chosen simply so that the sum of all filter coefficients 
is zero, in view of a similar condition on F (t) (equation 9). 



■13- 



The weights, w k , give some flexibility in the numerical approximation to the 
singular integral (equation 10). Some common choices are: 

1. Ramachandran & Lakshminarayanan (1971) 

w k = 2 for k odd, and w k = for k even (16) 

2. Shepp & Logan (1974) 

w k = 4k 2 /(4k 2 - 1) (17) 

3. Horn (1976) 

w k = 1 ( 18 ) 

The third set of weights corresponds to the trapezoidal rule for numerical 
integration or quadrature. Linear combinations of the above weights may also 
be used. For example, a combination of (1/3) of the first set and (2/3) 
of the third set produces weights which are alternately (2/3) and (4/3). 
This corresponds to Simpson's well-known rule for numerical quadrature. 
The second set of weights on the other hand corresponds to a numerical in- 
tegration formula which takes into account the singular nature of the in- 
tegral being approximated, as will be shown later. 

By summing the series indicated (equation 15b), one finds that 
(si) z ? = ir 2 /2, 4 and tt 2 /3 for the three sets of weights suggested. 
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CHOICE OF WEIGHTS. 

This analysis differs from the standard derivation of the weights w, . 
These coefficients are commonly obtained by inverse Fourier transformation 
from a filter response designed in the frequency domain. Their differences 
are usually discussed in a somewhat ad hoc fashion in terms of the need to 
low-pass filter the projection data in order to avoid aliasing or under- 
sampling. Clearly, this is wrong, since to avoid the effects of under-samp- 
ling, low-pass filtering has to be performed before sampling. After sampling 
we throw out the good with the bad, since they are no longer distinguishable. 
(Fortunately, the finite size of the detectors and to some extent the finite 
size of the source of radiation, account for some low-pass filtering of the 
projection data before sampling and thus help to limit the magnitude of the 
resulting artifacts). 

The derivation of these weights as coefficients in formulae for numeri- 
cal quadrature instead seems more insightful. The connection between these 
two points of view is made by Hamming [23, 24] in his discussion of the 
frequency response of integration formulae. 

Different choices of weights lead to different approximation behavior. 
As one might expect, there is a trade-off between noise and resolution. Ran- 
dom additive noise in the density integrals leads to noise in the final re- 
construction. The amplification factor depends on not only the details of 
scanning geometry (number of projections and number of rays per projection), 
but also the weights chosen. The first filter above (equation 16), for ex- 
ample, has fine resolution at the cost of sensitivity to noise and sharp 
contrasts, makes full use of the sampled projection and does not attenuate 
higher frequencies. The third filter, on the other hand (equation 18) lies 
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at the other extreme and tends to blur sharp edges, while suppressing noise; 
it removes some of the higher frequency components of the sampled projection 
data. The second filter (equation 17) lies between the two extremes. In 
practice, one should allow for the possibility of using different weights 
to suit different applications, in order to be able to exploit fully the 
trade-off between noise amplification and resolution. 

Overshoot in regions where p(£,e) does not have a continuous derivative 
with respect to £ is a common problem with filters that produce high resolu- 
tion results. They are most sensitive to violations of the assumptions un- 
derlying Radon's inversion formula. 

Finally, note that the two summations (equations 12 and 13) allow us to 
evaluate the estimated density at arbitrary points (r,$). In practice, one 
uses a fixed grid of picture cells, in the form of some regular tesselation 
of the plane. This limits the amount of computation and reflects the fact 
that resolution is limited, in any case, by the sampling width (si) along 
each projection, and that no new information is gained by performing re- 
construction on a grid much finer than this. 
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f~\ RAY-SAMPLING DENSITY. 



With the parallel-ray scheme described, sampling is uniform in i and e. 
That is, successive rays in a particular projection correspond to evenly 
spaced values of £, while successive projections correspond to evenly spaced 
values of e. Thus % and e are natural coordinates for the rays. Other co- 
ordinates are preferred when we are dealing with fan beams or more general 
scanning schemes. Essentially, whatever the scanning scheme, we must find 
coordinates £ and n natural to the particular geometry, such that we have 
uniform sampling in £ and n . The collection of ray sums pU,e) for a fixed 
value of n will be referred to as a generalized projection. It is simple now 
to rewrite the reconstruction formula (equation 8) as follows: 

1 l1m f f 
f(r,*) = — e + F (t) p(£,e) J d? dn (19) 

where, 

is the Jacobian of the transformation from (s,rt) space to (&,e) space. It 
can be conveniently visualized as the factor by which a small area in (5,n) 
space is expanded when mapped into U,e) space (see figure 5). 

Since we have uniform sampling in (?,n) space, the sampling density 
in (£,e) space equals the uniform density divided by J. To see this more 
clearly, let two rays (£,e) and U',e') be considered "near" each other if 
|* - A' | < &SL/2 and |e - e'| < 66/2. Clearly then the number of rays "near" 
a given ray is proportional to (1/J) 6£ 69. 
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%* Consequently, we can state that the ray-sampling 
density is inversely proportional to J. 

•3fr Further, it is clear that the contribution of a 
particular ray sum to a particular point must be 
weighted by J, that is, the inverse of the ray- 
sampling density. 

Intuitively, this seems reasonable since we do not want to emphasize 
contributions from regions of (£,0) space which happen to be sampled more 
densely than others. It should be noted that we can no longer expect all 
regions of the reconstruction to be equally well determined or resolved, 
since rays important to the reconstruction of one may be sampled more 
coarsely than the others. Fortunately, for practical fan-beam systems, 
/"■> the equivalent change in point-spread function over the region being recon- 
structed tends to be fairly small and thus not visually noticeable. 
We may write (equation 19): 



f(r,4>) = — / g(r,(j>,n) d n (21) 

4lT 2 J 



where 



lim f 

g(r,<j>, n ) = ^ J F e (t) J(?,n) p(g,n) d£ (22) 

In the general case, t = % - i l will be a function of both (r,(j>) and U,n). 
As a result, the inner integral may have to be evaluated separately for every 
point (r,<j>) in the reconstruction, for every projection. That is, g is a 
function of three variables, unless we further restrict the possible scanning 
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schemes. Fortunately, in most interesting cases a variable x(r,<fr) can be 
introduced which is natural to the scanning scheme such that g becomes a 
function of x and n only. (This may require splitting the variable t into 
a product of a term which depends on s and one which does not — the latter 
term can be moved out of the inner integral. 

If g can be written in terms of x and n only, a great computational 
efficiency arises, because the inner integral has to be evaluated only for 
every x(r,<t>) for a given projection, not separately for every picture cell 
(r,<|>). An example later on will make this clear. Frequently, a good choice 
for x(r,4>) is £' defined by the equation &(?')= V . 
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GENERAL LINEAR OPERATORS, 

If we can find a new parameter x(r,<j>) as described above, then the in- 
ner integral (equation 22) becomes 
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g(x.n) = /F £ [t(x,5,n)] J(5,n) p(e,n) d? (23) 

If we consider n as a parameter for the moment we can write this in a form 
that is more easily recognized: 

9 n (x) = y\(x.5) P n (5) & (24) 

This is a general linear operation with kernel K (x»?) = F (t) J. This 
operation is very similar to a convolution aside from the fact that in a 
convolution the kernel would be invariant. The above integral may also be 
referred to as a superposition integral and the general linear operator may 
also be called a linear space-variant operator. Integrals of similar form 
occur in the solution of partial differential equations, in which case the 
kernel is called a Green's function. In a number of special cases, such as 
uniform, parallel -ray scanning, the kernel is space-invariant (that is, is 
a function of x - £ only) and the operation simply becomes a convolution. 
Note, by the way, that the sampling-density factor, J, presents no 
special problems, representing merely a pre-multiplication of the ray sums. 
In fact, under fairly general conditions, J is a function of ? only and so 
each ray sum is simply multiplied by a factor depending on its position within 
its generalized projection. 
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PARALLEL-RAY SCANNING WITH VARIABLE RAY- AND PROJECTION-ANGLE SPACING. 

As an illustration of the utility of the new method for finding recon- 
struction algorithms, we develop an algorithm suited to parallel -ray scanning 
where both the spacing between successive rays in a projection and the in- 
terval between successive projection angles are non-uniform. 

Let the rays be evenly spaced in s, while projections are evenly spaced 
in n. Then we write £ as a function of § , and we write e as a function of 
n- Clearly, £(5) and e(n) should be monotonically increasing, continuous 
and differentiate. This also assures us that the inverse functions will 
exist. That is, given 1 we can find 5, and given e we can find n- The 
Jacobian (equation 20), here simplifies to: 

J=^iii (25) 

It is clear given these assumptions that J will be positive and that its 
two factors may be split between the inner and outer integrals. Now choose 
5', such that 2(5') = V , that is (equations 1 and 2), 

£(£') = r cos [e(n) - <f>] (26) 

Then, t = £(5 ) ~ £(5') and consequently we find that the inner integral 
is a function of 5' and n only. Finally (from equations 21 and 22), 

f(r,*) = -L- / g ( ? i, n ) |9. dn (27) 

4rr 2 J 9ri 



where 



lim f 
9(5 '' n) = e+oy F e <* - *'> Pfe»n) ff" <* (28) 
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Here then we have an inner integral which corresponds to a general linear 
operation. It becomes a convolution only if i happens to be a linear func- 
tion of £, that is, when the spacing of the rays in a given projection is 
uniform. 

For discrete sampling of the ray sums we approximate the above in- 
tegrals by sums: 



f(r,*) = — E g-U') 66. (29) 

47T 2 j J 3 



g 1'j" 2: F 1I1P1J 51, (30) 



Here again, p.. is the i ray sum in the j projection. If e. is the j 

f* k 

projection angle and a. is the distance of the i ray from the center of the 
region being scanned, then 66. is the angular interval associated with a 
particular projection, while Si. is the projection interval associated with 
a particular ray, where 

59 j = ( Vi " 9 j-i )/2 and 6£ i = (£ i+l " Vl )/2 < 31) 



Also, 



w . 
F = lU for i f V (32a) 

F.,.,64. , = E F... 81. (32b) 



That is, F . • is chosen so that 
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Here we happen to calculate g i , . for a set of values of V which corresponds 
to the set of values of % for rays whose ray sum is known. One could equally 
well have decided to perform the calculation for a different, perhaps evenly 
spaced set of values of A'. In either case, the values g.(£') have to be 
found from the known g i , by interpolation as indicated before. 

Note that the inner sum (equation 30) is not a convolution, but a general 
linear sum. Fortunately, it requires little more calculation than a simple 
convolution. It is also clear how the above simplifies if either the ray 
spacing or the projection-angle spacing becomes uniform. 

It should be pointed out that there is a minor practical problem due 
to the slow convergence of the series (equation 32b) for F.,.,. When rays 
are spaced evenly, this sum can be found analytically (equation 15b), while 
it is likely that numerical techniques are required here. If a is asymptotically 
linearly related to K, then the error term of the sum evaluated with n terms 
is proportional to 1/n. This illustrates the problem as well as suggests 
a solution. If we let s p be the sum of F...6A. from i = i'-n to i = i'+n, 
then a good estimate of the sum from i = -« to i = +~ is given by 

s. = (n+1) s n+1 - (n) s n (34) 
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TAKIN6 INTO ACCOUNT THE SINGULAR NATURE OF THE INNER INTEGRAL. 

So far, when we approximate the inner integral (equation 28) by the 
sum (equation 30), we pay little heed to the singular nature of the kernel 
FJa - £'). It is reasonable to suppose that better approximations may be 
found by considering methods which deal with the singularity. In this re- 
gard we note first that the values of the density integral p(5,n) are known 
only at discrete points, and that it is reasonable to assume that this com- 
ponent of the integrand varies relatively slowly over small distances along 
a projection (in fact, the partial derivative of this function was assumed 
to be continuous). The term (-1/t 2 ) on the other hand is known everywhere 
but varies rapidly near the point t = 0. We can make use of these observa- 
tions after splitting the inner integral into many integrals, each over the 
width of one detector. 

Let the i detector intercept rays lying in the beam between %. and 
i. + 1 (see figure 6). It measures the density integral p.. in the j th projec- 
tioru Further, let 

% { = U r + v +1 )/2 ( 35a ) 

and 



e = K+l " £ ., )/2 (35b) 



I * 

The inner integral (equation 28) can then be written as 



[Note that here the center of the i detector lies at (a. + Jt. +1 )/2], 



a.-.i £ i'+i IT IT. 



i'-l /i+1 1 p'+l +» *i+l 

Z J ("-) PU.e) dz + J [U p( £ , 9 ) di + £ J (- 1-) P U,e) dv 

^ i= - ro *i i , e i=i- + i . t* 
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Now note that t = (z - V) and 



/ 



£ i+l 



J- di- ] 



u i+1 - v 



£ i U"*') 2 U 1+1 - *') U t -*') (t 1+1 - £')U. - £') 



— (37) 



If we use this result and replace p(£,e) with p ., when a.s ^ < . then 
the inner integral (equation 36) becomes 



g i'j = 7: -7T p i'j 



(£,. fl - A,,) 






TJ 



(38) 
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The reader may verify that this in fact reduces to Shepp and Logan's method 
(equation 17) when rays are evenly spaced [that is, when i. = (i - 1/2)81 - w/2, 
6£ = W/(N - 1)]. More complicated integration formulae may be developed if 
one fits low-order polynomials to the values of p.. instead of assuming that 
the density integrals are constant over the width of one detector. Techniques 
for doing this may be found in standard texts on numerical analysis [23]. 
Here, however, we will be satisfied with the simple form developed above 
(equation 38) which is better than the method developed earlier 
(equation 30) since there is no difficulty in finding the weights by which 
the ray sums are to be multiplied. 

The method can be further simplified by using the other form (equation 
37) of the integral of (-1/t 2 ): 



i'-l 



i=-« L 



1 



1 



U i+1 -£') U.J-A') 



+0O 



p. . + 



( wy 



*y + 2 



i=i'+l 



1 



1 



U 1+1 -A') (*.-£') 



1J 



/*N 



(39) 
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Splitting both sums and rearranging terms leads to the surprisingly simple 
result, 

The reader may want to compare this with the original form of the inner in- 
tegral (equation 4), from which this result can also be obtained directly. 
This numerical approximation of the inner integral is particularly advanta- 
geous from a computational point of view since it is no longer necessary to 
keep a two-dimensional array of pre-calculated weights. This assumes that 
one can afford to calculate {i. -a'), and that the ray sums are replaced 
by the differences of ray sums as required above (equation 40). This latter 
calculation is needed only once per projection. 

The form of the result also implies that reconstruction may be possible 
when the ray spacing varies discontinuously, that is, when i(z) does not 
have a continuous derivative with respect to £. This in turn suggests the 
possibility of using evenly spaced detectors; combining the measurements of 
neighboring detectors in portions of the projection where high resolution 
is not required. 
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MOTIVATION FOR STUDYING VARIABLE RESOLUTION METHODS. 
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In a number of situations, one is intested in reconstructing as object 
buried inside some larger entity of less interest. If one were simply to 
restrict the scanned region to the object of interest, correct reconstructions 
would not be obtained, since the absorbing density is then non-zero outside 
the region being scanned. This violates assumptions underlying Radon's 
formula. Up to now, the only alternative was to scan the whole region oc- 
cupied by absorbing material and reconstruct it with uniform resolution. 
(At best, there is some saving in the backprojection step, since one need 
not calculate the density of picture cells outside the region of interest). 

The new variable resolution method to be illustrated here has the ad- 
vantage that the computation of the filtered projection is speeded up con- 
siderably since fewer ray sums have to be measured. Of equal importance may 
be the fact that less radiation is needed to obtain this smaller set of 
measurements. 
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DEMONSTRATION OF THE VARIABLE RESOLUTION METHOD, 

In order to illustrate some of the features of the new method, a computer 
algorithm based on the results derived here (equations 40 and 29) has been 
developed. This algorithm has been applied to ray sums calculated from a 
mathematically defined object or phantom composed of elliptical parts (see 
figure 7). A comparison will be made of the results obtained in three 
cases: 

(a) .N =200 evenly spaced rays, 2 mm apart. 

(b) N = 100 evenly spaced rays, 2 mm apart at the center, 
8 mm at the edge. 

(c) N = 100 evenly spaced rays, 4 mm apart. 

In order that the comparison be fair, all other parameters were held constant. 
The region scanned had a diameter W = 400 mm, and M = 150 projection angles 
were employed in each case. For case B, the following transformation from 
uniform scanning coordinates to actual scanning coordinates was used: 

l =|(f) (1 + S 2 ) (41) 

where -1 - £ ^ +1 . In all cases the ray sums were averages obtained by in- 
tegrating from the left edge of a beam striking a detector to the right edge 
of this beam so as to simulate the suppression of high frequency components 
obtained in practice as a result of the finite width of the detectors. 

Each projection is first processed to produce the differences required 
in the summation (equation 40). The filtered sums are then determined for 
positions corresponding to the individual detectors. In order to facilitate 
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back-projection, these values are (linearly) interpolated to a much finer, 
evenly spaced set of positions. Backprojection proceeds much as it does 
for convolutional algorithms by processing each picture cell in turn. For 
every picture cell, the appropriate point in the interpolated filtered data 
is found by considering the projection angle (as illustrated in figure 4). 
The value found there is then added into the sum accumulated for this pic- 
ture cell so far. The whole process is repeated for all projection angles. 

The picture cells were spaced 1.5 mm and lay inside a circle of diameter 
330 mm for the reconstructions shown in figures 8 and 9. For the reconstruc- 
tions shown in figures 10 and 11, the spacing was .75 mm inside a circle of 
150 mm diameter. 

The time required for backprojection was essentially the same for the 
three cases, while the time for the general linear operation in cases B and 
C was about one quarter that required for case A, as expected. 
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ANALYSIS OF RECONSTRUCTION RESULTS. 



The low resolution near the edge of the region of reconstruction (where 
the rays are spaced 8 mm apart) of method B can best be seen in figure 8. 
Hot much is visible in this figure of the central components however. Figure 
9 more clearly shows the low overall resolution of method C. (It is important 
not to be misled by the apparent high resolution of high contrast features 
due only to the reduced density scale of this mode of presentation). The 
good resolution of method B in the central regions is illustrated by figure 
10, as well as by the density profiles in figure 11. The density profiles 
are along the lines indicated in figure 12. It appears that while the 
variable resolution method, B, requires only about as much computation as 
the low resolution method, C, it has about as much central resolution as 
the high resolution method, A. 

The reader will have noticed the reconstruction artifacts particularly 
apparent in the high contrast presentations (figures 9 and 10). The phantom 
was purposefully constructed to include high contrast features outside a 
central region with a variety of low contrast features, since artifacts 
radiating outwards from the former often degrade the presentation of the 
latter. 

As indicated earlier, these artifacts are due to the projection data's 
failure to obey the assumptions underlying Radon's formula. That is, the 
partial derivative of pU,e) with respect to i is not everywhere continuous. 
It is easy to see that the discontinuities occur at the edges of the elliptical 
components of the phantom and that line-like artifacts oriented tangentially 
to the high-contrast components radiate across the reconstructed density 
distribution. The exact magnitude of these artifacts depends on the particular 
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alignment of a projected edge relative to the edges of the detectors. It is 
easy to see, too, that the magnitude of these effects varies inversely with 
the number of projection angles, since the contribution of each projection 
varies in this way. Further, it can be shown that the magnitude of this ef- 
fect also decreases with the number of rays in a projection. The artifacts 
are visible in the examples presented here because both the number of pro- 
jection angles (150) and the number of rays per projection (100 or 200) are 
relatively small and because of the strong contrast between some of the large 
features in the phantom. 

The important point is that these artifacts, while visible, do not 
mask the low contrast detail in the center, and that the magnitude of the 
artifacts in the central region are not significantly larger in the recon- 
struction obtained using method B, than they are in those obtained using 
method A. The new method would be of less interest if this were not the 
case. 



/""N 



-31- 



SUMMARY AND CONCLUSIONS. 

A technique has been developed for finding reconstruction methods for 
arbitrary ray-sampling schemes. The algorithms use a general linear opera- 
tor, the kernel of which depends on the particular scanning geometry. It 
is suggested here that the kernel coefficients can be fruitfully considered 
as weights in a method for numerical quadrature of an integral. In a few 
special cases, the kernel is a function of the difference of coordinates 
only, and the general linear operation becomes simply a convolution. In 
this case, the algorithm essentially reduces to the familiar convolution- 
backprojection-summation method. 

As an illustration of the more general case, an algorithm was derived 
which applies to parallel -ray data when the spacing of rays in a projection 
is uneven and the projections are spaced unevenly in angle. The new method 
requires little more computation than does the convolutional method. 

Importantly, the Fourier transform is not used in the derivation. In- 
deed, it does not apply to the linear space-variant systems that occur in 
the general case, and also in particular cases of practical importance. 
A future paper will explore the application of this general method to a 
variety of fan-beam scanning geometries suitable for modern tomographic 
machines. 
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FIGURE CAPTIONS. 

Figure 1. Parallel-ray scanning geometry. Many projections are measured, 
each with rays arriving in a particular direction. 

Figure 2. Fan-beam scanning geometry. Many projections are measured, each 
with the source in a particular position. 

Figure 3. Designation of particular rays and calculation of distance be- 
tween a ray and a point in the region being scanned. 

Figure 4. Detailed geometry of a parallel ray projection. 

Figure 5. Transformation from uniform scanning coordinates to coordinates 
used in Radon's inversion formula. The Jacobian is the ratio of 
the area of the quadrilateral A'B'C'D 1 to that of ABCD. 

Figure 6. Positions of detectors along projection as defined in derivation 
of numerical approximation to singular integral. 

Figure 7. Outlines of the elliptical components of the phantom used in the 
demonstration of the variable resolution method. The numbers in- 
dicate the absorbing densities of the components. 

Figure 8. Reconstructions obtained in the following three cases: 

(A) 200 evenly spaced rays 

(B) 100 unevenly spaced rays 

(C) 100 evenly spaced rays 

In this figure black corresponds to a density of -.06, while white 
corresponds to a density of 1.22. 
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Figure,9*. The same reconstructions as in previous figure displayed with higher 
contrast. Here black corresponds to a density of .94, while white 
corresponds to a density of 1.06. 

Figure 10. Higher resolution display of central regions of the reconstructions 
shown in previous figures. 

Figure 11. Horizontal density profiles through the middle of the central 

circular component of each of the three reconstructions. While 
method B requires only about as much computation as method C, it 
gives rise to resolution about as good as method A. 

Figure 12. Lines along which density profiles of previous figure were taken. 
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